clear all
close all
clc

global Glob 
%%% Monte Carlo SIMULATIONS


% Size of markets
dimNum=30;
dimNum2=20;
dimS=600;
% True parameters
sigma0=0.25;
omega0=[0.85;0.9;0.95];
dimProv=length(omega0);
nationGrid=[1;2];

dimNation=length(nationGrid);
alpha0_1=[1 0.97 ;0.9 1];
alpha0_2=[1 0.9 ;0.97 1];
alpha0_3=[1 0.8 ;0.85 1];
alpha0=cat(3,alpha0_1,alpha0_2,alpha0_3);
seed0=20212025;
dimMC=250;
maxiter=500;

param.alpha=alpha0;
param.omega=omega0;
param.sigma=sigma0;
Glob.dimS=dimS;
Glob.dimProv=dimProv;
Glob.dimNum=dimNum;
Glob.dimNation=dimNation;
Glob.seed=seed0;

ind=logical(1-eye(dimNation));
ind=cat(3,ind,ind,ind);
theta0=[omega0;sigma0;alpha0(ind)];

Glob.dimNum=dimNum2;
tic
TrueMarriage=solveMonteCarlo20191030(param,Glob)
toc
disp('Number of Marriages:')
q=TrueMarriage;sum(reshape(q,2*2,3))


Glob.TrueMarriage=TrueMarriage;
Glob.dimNum=dimNum;

STOREPARAM=zeros(length(theta0),dimMC);
STOREFCT=zeros(dimMC,1);

for iterMC=1:dimMC
disp(['Iteration: ' num2str(iterMC)] )

%% Solving the true allocation
c=clock;
myseed=ceil(c(6)*10000);
disp(myseed)
Glob.seed=myseed;

tic
fvalM0=estimMonteCarlo20191030(theta0);
toc

% Estimation 
options=optimset('Display','none','MaxIter',maxiter,'MaxFunEvals',10000,'TolFun',1e-6,'TolX',1e-4);
step=0.05;
[xM,fvalM1,exitflag,output,vinit,fvinit]=fminsearch_j(@estimMonteCarlo20191030,theta0,options,step,0.1);                
disp(fvalM1)

[xM theta0]
STOREPARAM(:,iterMC)=xM;
STOREFCT(iterMC)=myseed;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dimGrid=100;

F=zeros(dimGrid,length(theta0));
XI=zeros(dimGrid,length(theta0));
MIN=theta0'-max([theta0'-min(STOREPARAM');max(STOREPARAM')-theta0']);
MAX=theta0'+max([theta0'-min(STOREPARAM');max(STOREPARAM')-theta0']);
MIN=floor(MIN*100)/100;
MAX=ceil(MAX*100)/100;
for i=1:length(theta0)
    stepxi=(MAX(i)-MIN(i))/(dimGrid-1);
    xi=MIN(i):stepxi:MAX(i);
    [f,xi] = ksdensity(STOREPARAM(i,:)',xi');
    F(:,i)=f'/sum(f);
    XI(:,i)=xi';
end
STOREPARAM1=STOREPARAM;

F2=zeros(dimGrid,length(theta0));
for i=1:length(theta0)
    stepxi=(MAX(i)-MIN(i))/(dimGrid-1);
    xi=MIN(i):stepxi:MAX(i);
    [f,xi] = ksdensity(STOREPARAM(i,:)',xi');
    F2(:,i)=f'/sum(f);
end

STOREPARAM2=STOREPARAM;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FIGURE F1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=figure('Name','FIGURE F1');
tiledlayout(3,4);
p1=nexttile; % omega_1
plot(XI(:,1),F(:,1),'LineWidth',2,'Parent',p1)
hold on
plot(XI(:,1),F2(:,1),'LineWidth',2,'Parent',p1)
hold on
y = ylim; % current y-axis limits
plot([omega0(1) omega0(1)],[y(1) y(2)],':','LineWidth',3,'Parent',p1)
title('$\omega_A$','FontSize',24,'interpreter','latex')
ylim(p1,[y(1) y(2)])

p2=nexttile; % omega_2
plot(XI(:,2),F(:,2),'LineWidth',2)
hold on
plot(XI(:,2),F2(:,2),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([omega0(2) omega0(2)],[y(1) y(2)],':','LineWidth',3,'Parent',p2)
title('$\omega_B$','FontSize',24,'interpreter','latex')
ylim(p2,[y(1) y(2)])

p3=nexttile; % omega_3
plot(XI(:,3),F(:,3),'LineWidth',2)
hold on
plot(XI(:,3),F2(:,3),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([omega0(3) omega0(3)],[y(1) y(2)],':','LineWidth',3,'Parent',p3)
title('$\omega_B$','FontSize',24,'interpreter','latex')
ylim(p3,[y(1) y(2)])

p4=nexttile; % sigma_1
plot(XI(:,4),F(:,4),'LineWidth',2)
hold on
plot(XI(:,4),F2(:,4),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(4) theta0(4)],[y(1) y(2)],':','LineWidth',3)
title('$\sigma$','FontSize',24,'interpreter','latex')
ylim(p4,[y(1) y(2)])

p5=nexttile; % alpha_12^1
plot(XI(:,5),F(:,5),'LineWidth',2)
hold on
plot(XI(:,5),F2(:,5),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(5) theta0(5)],[y(1) y(2)],':','LineWidth',3)
title('$\alpha_{12}^A$','FontSize',24,'interpreter','latex')
ylim(p5,[y(1) y(2)])

p7=nexttile; % alpha_12^2
plot(XI(:,7),F(:,7),'LineWidth',2)
hold on
plot(XI(:,7),F2(:,7),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(7) theta0(7)],[y(1) y(2)],':','LineWidth',3)
title('$\alpha_{12}^B$','FontSize',24,'interpreter','latex')
ylim(p7,[y(1) y(2)])

p9=nexttile; % alpha_12^3
plot(XI(:,9),F(:,9),'LineWidth',2)
hold on
plot(XI(:,9),F2(:,9),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(9) theta0(9)],[y(1) y(2)],':','LineWidth',3)
title('$\alpha_{12}^C$','FontSize',24,'interpreter','latex')
ylim(p9,[y(1) y(2)])

p6=nexttile(9); % alpha_21^1
plot(XI(:,6),F(:,6),'LineWidth',2)
hold on
plot(XI(:,6),F2(:,6),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(6) theta0(6)],[y(1) y(2)],':','LineWidth',3)
title('$\alpha_{21}^A$','FontSize',24,'interpreter','latex')
ylim(p6,[y(1) y(2)])

p8=nexttile(10); % alpha_21^2
plot(XI(:,8),F(:,8),'LineWidth',2)
hold on
plot(XI(:,8),F2(:,8),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(8) theta0(8)],[y(1) y(2)],':','LineWidth',3)
title('$\alpha_{21}^B$','FontSize',24,'interpreter','latex')
ylim(p8,[y(1) y(2)])

p10=nexttile(11); % alpha_21^3
plot(XI(:,10),F(:,10),'LineWidth',2)
hold on
plot(XI(:,10),F2(:,10),'LineWidth',2)
hold on
y = ylim; % current y-axis limits
plot([theta0(10) theta0(10)],[y(1) y(2)],':','LineWidth',3)
title('$\alpha_{21}^C$','FontSize',24,'interpreter','latex')
ylim(p10,[y(1) y(2)])

l=legend({'2000','20000'},'FontSize',18);
leg_pos = get(l,'position') ;
title(l,{'Number of','Simulations'},'FontSize',18)
set(l,'position',[0.73,0.25,0.25,0.2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% TABLE F1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
name=['$\omega_A$     ';'$\omega_B$     ';'$\omega_C$     ';'$\sigma$       ';'$\alpha_{12}^A$';'$\alpha_{21}^A$';'$\alpha_{12}^B$';'$\alpha_{21}^B$';'$\alpha_{12}^C$';'$\alpha_{21}^C$'];
name=['$\omega_A$       ';'$\omega_B$       ';'$\omega_C$       ';'$\sigma$         ';'$\alpha_{A}(I,F)$';'$\alpha_{A}(F,I)$';'$\alpha_{B}(I,F)$';'$\alpha_{B}(F,I)$';'$\alpha_{C}(I,F)$';'$\alpha_{C}(F,I)$'];
Sep=repmat(' & ',[length(theta0) 1]); 
End=repmat('\% \\ ',[length(theta0) 1]);
SepP=repmat(' \% & ',[length(theta0) 1]); 
disp([name Sep num2str(theta0,2) Sep Sep num2str((mean(STOREPARAM1,2)-theta0)./theta0*100,3) SepP num2str(std(STOREPARAM1,[],2)./theta0*100,3) SepP Sep  num2str((mean(STOREPARAM2,2)-theta0)./theta0*100,3) SepP num2str(std(STOREPARAM2,[],2)./theta0*100,3) End])










